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We investigate the onset of collective oscillations in a network of pulse-coupled leaky-integrate- 
and-fire neurons in the presence of quenched and annealed disorder. We find that the disorder 
induces a weak form of chaos that is analogous to that arising in the Kuramoto model for a finite 
number N of oscillators [O.V. Popovich at al, Phys. Rev. E 71 065201(R) (2005)}. In fact, the 
maximum Lyapunov exponent turns out to scale to zero for N — > oo, with an exponent that is 
different for the two types of disorder. In the thermodynamic limit, the random-network dynamics 
reduces to that of a fully homogenous system with a suitably scaled coupling strength. Moreover, 
we show that the Lyapunov spectrum of the periodically collective state scales to zero as l/N 2 , 
analogously to the scaling found for the 'splay state'. 
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I. INTRODUCTION 



One of the most general and relevant dynamical phenomena observed in the mammalian brain is the rythmic 
coherent behaviour involving different neuronal populations The dynamics of neural circuits has been widely 
studied, by invoking various kinds of neuron models; collective oscillations are commonly associated with the inhibitory 
role of interneurons 0. However, coherent activity patterns have been observed also in "in vivo" measurements of 
the developing rodent neocortex and hyppocampus for a short period after birth Q , despite the fact that at this early 
stage the nature of the involved synapses is essentially excitatory 0] • 

Independently, theoretical studies of fully coupled excitatory networks of leaky integrate-and-fire (LIF) neurons have 
revealed the onset of macroscopic collective periodic oscillations (CPOs). This dynamical state is quite peculiar: 
the collective oscillations are a manifestation of a partial synchronization among the neuron dynamics and this is 
one way of identifying this phenomenon, which is, however, more subtle: the macroscopic period of the oscillations 
| does not coincide with (is longer than) the average interspike-interval ISI of the single neurons and the two quantities 
are irrationally related. In fact, this phenomenon is also called self-organized quasi periodicity and can be observed 
in a wide class of globally coupled systems Q- In the context of pulse-coupled neural networks, CPOs arise from 
the destabilization of a regime characterized by a constant mean-field and a strictly periodic evolution of the single 
' neurons: this regime, termed "splay state", has been widely studied in several contexts, including computational 
^\ [ neuroscience Q. 

0^ ' Since real neural circuits are not expected to have a full connectivity [9|], it is important to investigate the role of 
dilution on the occurrence of the stability of CPO. We do so by investigating an excitatory network of LIF neurons 
^1 . with 20% of missing links in two different setups: quenched disorder, where the network topology is fixed; annealed 
disorder, where the active connections are randomly and independently chosen at each pulse emission. As a first step, 
we rewrite the dynamical equations as a suitable event driven map, by extending the approach developed in [l(| • We 
do so by introducing a pair of variables for each neuron, to account for the evolution of the local electric field. This 
step is particularly important for the computation of the Lyapunov exponents, as it allows expressing the evolution 
equations into a "canonical" form and thereby simplifies the implementation of standard dynamical-system tools. 

We find that the regime of CPOs is robust against the presence of dilution, both in the quenched and annealed 
setup. However, at variance with the homogeneous fully-coupled case, the dynamics of finite disordered networks 
turns out to be chaotic, although the degree of chaoticity decreases with the number N of neurons. In fact, the 
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maximum Lyapunov exponent goes to zero as 1/iV 3 . The exponent /3 is smaller in the quenched setup, indicating 
that hnite-size effects are stronger. In the homogeneous case, we are able to determine the full Lyapunov spectrum 
for sufficiently large numbers of neurons. As a result, we find that the first band of the spectrum scales as that of the 
splay state [a EL namely as 1/N 2 

The paper is organized as follows. In section [TTI we introduce the model and the event-driven map that is used to 
carry out the stability analysis. In Sec. IIIII we discuss the collective dynamical behaviours observed in the presence 
of disorder. Sec. IIVI is devoted to the Lyapunov analysis of these coherent solutions in the large N limit. Finally, in 
Sec. El we summarize the main results and the open problems. 



II. THE MODEL 

We study a network of N neurons, whose individual dynamics is modelled as a LIF oscillator. Following Refs. [To| . 
the membrane potential Xi (t) of the i — th oscillator evolves according to the differential equation 

x i (t) = a-x i (t)+gE i (t) i = l,---,N (1) 

where a > 1 is the suprathreshold input current and g > gauges the coupling strength of the excitatory interaction 
with the neural field Ei (t) . At variance with the fully-coupled network, where all neurons depend on the same "mean 
field" E(t), here we consider a general setup, where neurons have different connectivities inside the network. As 
a result, it is necessary and sufficient to introduce an explicit dependence of the neural field on the index i. The 
discharge mechanism operating in real neurons is modelled by assuming that when the membrane potential reaches 
the threshold value X4 = 1, it is reset to the value Xi = 0, while a pulse is transmitted to and instantaneously received 
by the connected neurons. The field Ei(t) is represented as the linear superposition of the pulses received by neuron 
i at all times t n < t: the integer index n orders the sequence of the received pulses. Each pulse is weighted according 
to the strength of the connection Cj^i between the emitting (j(n)) and the receiving (i) neuron. In general, the 
connectivity matrix C is assumed to be non-symmetric. In formulae, 

^(*)= E C 3(n)At ~tn)e(t- t<A ex P [-«(*-*„)] (2) 

n\t n <t ^ ' 

where 6{x) is the Heavyside function and a is the inverse pulsewidth (the pulse shape has been chosen as in Ref. . 
Since also in the diluted case we will consider massively connected networks [121 ] , where the average number of synaptic 
inputs per neuron varies proportionally to the system size, it is natural to scale the synaptic strength with N as done 
in Eq. ©. 

The model is fully characterized by the two sets of equations ((T|) and |(5J). The dynamical system has a peculiar 
mathematical structure, with the field variable appearing as a memory kernel, which involves a summation over all 
past spiking events. We find much more convenient to turn the explicit equation @ into the implicit differential 
equation 

2 

Ei(t) + 2aEi(t) + o?E^t) =1 J C j(n)<i S(t - t n ) . (3) 

n\t n <t 

As a result, the dynamics of the neural network model takes the more "canonical" form of a set of coupled ordinary 
differential equations (|1|3|) . which can be analyzed with the standard methods of dynamical systems (see e.g., 5, 8. 

MM)- 



A. Event-driven map 



The presence of <5-like pulses into the set of coupled differential equations (fTJ) and © may still appear as an intrinsic 
technical difficulty for the estimation of the stability properties. Actually, the standard algorithms for the evaluation of 
Lyapunov exponents rely upon the integration of differentiable operators acting in tangent space (see [Hj]). However, 
one can easily get rid of this problem by tansforming the differential equations into a discrete time event-driven 
mapping. This task can be accomplished by integrating Eq. (J3J from time t n to time i„+i (t n representing the time 
immediately after the emission of the n-th pulse). An explicit mapping can be written by introducing the auxiliary 
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variable Qi :— aEi + Ei, 

Ei(n + 1) = E t {n)c- aT{ - n *> + Q 4 (n)r(n)e- QT(n) (4a) 

2 

Qi{n + 1) = Qi(n)e- aT W + C j(n+1) ~ (4b) 

Xi(n + 1) = x l (n)e^ T(n) + <z(l - e - r(Tl) ) + gH^n) . (4c) 

Here r(n) = t n+ \ — t n is the interspike time interval: it is determined by the largest membrane potential (identified 
by the label m(n)) reaching the threshold value x m — 1 at time i„+i, 



r(n) = In 

where 



a - x m (n) 



a + gH m (n) - 1 



(5) 



a— 1 \ a — 1 y (a — 1) 

for the parameter values considered in this paper {g > and a < 1), Hi{n) > 0. 

Altogether, the model now reads as a discrete time map of 3N — 1 variables, {Ei, Qi, x{\: one degree of freedom, 
x m (n) = 1, is eliminated as a result of having constructed the discrete-time dynamics with reference to a suitable 
Poincare-section. At variance with the usual approach (see e.g. Ref. 14]), the evolution equation does neither involve 
(5-like temporal discontinuities, nor formally infinite sequences of past events: the map is a piecewise smooth dynamical 
system. 

For the sake of simplicity, we assume that the entries of the connectivity matrix Cjj are either or 1 (the homoge- 
neous fully-coupled case correspondingto Cj t i = 1 for all j's and i's). If the entries are chosen randomly, symmetries 
are lost and the only way to determine the asymptotic state for a finite N is by numerically simulating map (U). In 
what follows, we consider two different setups: the synaptic connections are randomly chosen and are constant in 
time (quenched disorder); each time a neuron fires, the neurons receiving the excitatory pulse are randomly chosen 
(annealed disorder). 

B. Linear Stability Analysis 

As usual, the stability of (U) can be analyzed by following the evolution of infinitesimal perturbations in the tangent 
space. The corresponding equations are obtained by linearizing Q as follows, 

5Ei(n + 1) = e- aT ^5Ei{n) + T(n) e - aT ^6Q t (n) - (aE t (n)(aT(n) - l)Q l {n)) e~ aT{n) ST(n) , (7a) 
5Qi(n + 1) = e- aT ^ [SQi(n) - aQi(n)} 5r(n) , (7b) 
5xi(n + l) = er T{n \8xi{n) - Xi{n)5T{n)]+ ae^ T 5r{n) + g5Hi{n) i = l,...,N -1 ; 6x m {n + 1) = . (7c) 

An explicit expression of <5r(n) can be obtained by differentiating Eqs. (I5I6D 

5t(ti) = T x Sxi(n) + T E 5E(n) + tqSQ(ti) , (8) 

where t x := dr/dxi and analogous definitions are adopted for te and tq. Moreover, SHi(n) is a short-cut notation 
for the linearizazion of expression ©, which in turn depends on 5Ei{n),8Qi{n) and Sr(n). For more mathematical 
details see Q3 GUGJ]. 

The degree of chaoticity of a given dynamical state is obtained by computing the Lyapunov spectrum, i.e. the set 
of 3N — 1 exponential growth rates Xi along the independent directions in tangent space. The Lyapunov spectrum 
has been numerically estimated by implementing the standard algorithm (T3 |. 

III. COLLECTIVE DYNAMICS 



In the fully coupled homogeneous case the local fields Ei, Qi are independent of the index i and the number of 
equations reduces to N + 1. Depending whether a is smaller or larger than a critical value a c (g), the dynamics either 
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FIG. 1: Minimal and maximal values of E(n) for different a values for N = 1, 600 with quenched disorder and p = 0.4. 



converges to a so-called splay state, with constant E(t 
evolve periodically and quasi-periodically, respectively 

This "mean-field" dynamics is expected to change in disordered networks 
fields Ei and Qi, we find it convenient to introduce the average variables, 



or to a partially synchronized state, where E(t) and Xi(t) 

a. 

Given the neuron dependence of the 



1 N 



k=l 



1 * 



(9) 



fe=i 



as a tool to characterize the resulting dynamical regimes and to compare with the homogeneous case. All the results 
reported in this paper have been obtained by fixing the fraction / of missing connections equal to / = 0.2. Numerical 
investigation for different choices of / (not reported here) show similar results. In analogy with the analysis of the 
homogeneous case, we have chosen a as the main control parameter (one can easily realize that choosing g would be 
equivalent). 

In Fig. [1] we plot the maximum and the minimum values of E(n) for different values of a, for g — 0.4 and ./V = 1600 
in the presence of quenched disorder. The bifurcation diagram is similar to the one observed in the globally coupled 
networks [f|. However, the splay state found for small a- values has been replaced by a fluctuating asynchronous 
state, where the average field E(t) is only approximately constant (the difference between mm E(n) and maxE(n) is 
of the size of the symbols). The periodic collective state is analogously affected by small irregular fluctuations. This 
strong similarity between globally coupled and the diluted networks is not surprising: for any finite value of /, upon 
increasing N, the differences among the fields should progressively disappear. In fact, in the limit N — > oo, 

we expect randomly diluted networks to behave as fully coupled ones, provided that variables and parameters are 
properly rescaled. More precisely, the average field E(t) in a network with a fraction / of missing links and a coupling 
constant g, is expected to be equivalent to the neural field E generated in a fully coupled network with coupling 
constant g(l — /). This expectation is confirmed by the critical value a c separating the two dynamical phases. As 
shown in Fig. [TJ for g = 0.4 and / = 0.2, a c ~ 6.8, a value which coincides, within the statistical error, with the 
critical value found in a globally coupled network with g — 0.32 = 0.4 x (1 — 0.2) [ltj . A further consequence of the fact 
that the dynamics of disordered networks reduces in the thermodynamic limit to that of homogenous fully coupled 
ones is that the sample-to-sample fluctuations typical of quenched disorder should vanish as well. For this reason we 
have not performed averages of different realizations of the disorder. 

A more detailed representation of the quenched dynamics above a c is obtained by looking at the projection in the 
plane (E(n), Q(n)). In Fig. [2] data sets are shown for a = 9 and increasing values of N: panels a and b correspond 
to the annealed and quenched case, respectively. This allows seeing that the two kinds of disorder yield indeed 
qualitatively similar results. 

More precisely, the phase points cluster around closed curves, revealing a "noisy" periodic dynamics. In the annealed 
case, the fluctuations can be attributed to the stochasticity of the evolution rule. Surprizingly, they are even larger 
in the quenched case, which correspond to a deterministic dynamics, in which case they must be attributed to the 
presence of deterministic chaos. 

Upon increasing N, the amplitude of the fluctuations decreases and the attractor shape converges to an asymptotic 
curve, which should correspond to that of a homogeneous fully coupled network with a properly rescaled coupling 




FIG. 2: (Color online) E(n) versus Q(n) for various system syze: (a) annealed disorder and (b) quenched case here is shown 
also the annealed result for N = 100, 000. All data refers to a — 9. 



strength g. This expectation is indeed confirmed in Fig. [3J where the attractor of a homogeneous network (with 
N = 1600 and g — 0.32) superposes to that of the annealed dynamics (with N = 10 5 and g = 0.4). In the quenched 
case, the asymptotic shape is the same, but the convergence is slower. 




FIG. 3: gE(n) versus gQ(n). The (blue) curve represents the attractor of a globally coupled system with N = 1,600 and 
g = 0.32, while the black dots refer to the annealed case with N = 100, 000 and g = 0.4. For both cases a — 9 and a = 1.3. 



In order to investigate quantitatively the scaling behavior of the finite-size corrections, we have studied the N 
dependence of 

AQ = (Q)(N) - (Q)(oo) (10) 

where the angular brackets denote the (time) average of Q-value of all configurations falling within the lower gE(n)- 
window [0.36, 0.44]. Since the asymptotic value (Q)(oo) is independent of the setup, we have extrapolated it in the 
simpler context of a fully coupled network with g — 0.32 As a result, we find that AQ always converges to zero as 
power-law, N~° , with (3 = 1 in the fully coupled network, f3 = 0.55 for quenched disorder, and f3 — 0.82 for annealed 
disorder (see FigHk). These latter values have to be considered as approximate estimates, affected both by statistical 
errors and finite-size corrections. A more accurate estimate is beyond our computational capabilities. However, for 
the present purpose, the relative differences are large enough to conclude that quenched disorder is characterized by 
a slower convergence. 

As a further test of the collective dynamics, we have computed the standard deviation cr(n) of the instantaneous 
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FIG. 4: (Color online) (a) The quantity AQ (see Eq. I10|) versus the network size; (black) triangles refer to the quenched case, 
while (red) squares to the annealed one: data fits are reported as dashed lines. The reported values have been obtained by 
averaging over 6 x 10 6 to 3 x 10 7 consecutive spikes, (b) Standard deviation a (see Eq. Ill|) vs N : circles and diamonds refer 
to the annealed and quenched case, respectively. The standard deviation have been averaged over M — 900 iterates. All data 
refers to / = 0.2, a = 9, g = 0.4 and a = 1.3. 



fields 

^) = (^M-^») 1/2 , (id 

which is a measure of the degree of synchronization among the various fields. In Fig. HJd we plot the time average (a) 
for annealed and quenched disorder: in both cases (a) decreases with N as iV -1 / 2 . These results confirm once more 
that in the limit N — > oo the neural field dynamics converges to that of homogeneous networks, irrespectively of the 
disorder. 



IV. LYAPUNOV ANALYSIS 



In order to provide a more detailed characterization of the macroscopic as well as of the microscopic dynamics, in 
this section we analyse the Lyapunov spectra for the fully-coupled network and for its disordered variants. 



A. Globally Coupled Network 



In this case, the fields seen by the neurons are equal to one another and it is therefore sufficient to introduce a 
single pair of field variables E and Q. The corresponding stability analysis for the splay state has been analytically 
carried out in Ref. 10], finding that the spectrum of Floquet exponents is composed of a band of values of order 1/N 2 , 
which become of order 0(1) at one band extremum, plus two isolated eigenvalues associated with the field dynamics. 
Therefore, it is convenient to start the numerical analysis by testing the effect of attaching a pair of field variables 
to each neuron. The results plotted in Fig. [S£i refer to the Lyapunov spectrum of the splay state for a — 3 and 
g = 0.4: one can observe two bands and two isolated exponents. The first band, composed of N — 1 nearly vanishing 
exponents, and the two isolated exponents coincide with the eigenvalues analytically obtained in (Iol |. The second 
band, composed of 2(N — 1) exponents quantifies the transversal stability of the synchronization manifold Ei = E%, 
Qi = Qi Vi, where, without loss of generality, we consider the field of the 1st neuron as the reference variable. In 
fact, if one introduces the suitable difference variables, 

w i = E i -E 1 z i = P i -P 1 i = 2, • • • , JV . (12) 

the corresponding tangent dynamics is ruled by the equations 

5wi(n + l) = 8w t {n)e- aT ^ + 8z t {n)T{n)e- ar ^ (13) 
Szi(n + l) = 5z t {n)e- aT{n) i = 2,---,N . (14) 
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FIG. 5: (Color online) Lyapunov exponents Aft versus k/N . Complete spectrum for TV = 50 for a = 3 (a) and a = 9 (b). 
Rescaled Lyapunov eigenvalues of the first band Aft x TV 2 for different system sizes, namely TV = 50 (filled black circles), f00 
(open red squares), 200 (green open triangles) and 400 (blue crosses) , for a = 3 (c) and a — 9 (d). For a = 3 in (c) also the 
analytical expression (dashed black line) reported in Eq. (29) in [ll[ is shown. The data has been obtained by following the 

evolution in the tangent space of the event driven map for a number of consecutive spikes of the order of 10 8 10 9 , after 

discarding a transient composed by 50, 000 x TV spikes. 



This shows that the linearized dynamics of each pair of twin variables (5wi,Szi) is decoupled from the rest of the 
network. Its stability can be evaluated by solving the corresponding two-dimensional eigenvalue problem. Accordingly, 
we expect two bands of N — 1 eigenvalues each. However, since the two eigenvalues of Eq. (IT5|) are both equal to —a, 
we do have a single band, as indeed observed in Fig. [5^,. As a last check, we have verified the 1/N 2 scaling of the first 
band of the splay state. The nice overlap of the three sets of exponents corresponding to N = 100, 200, and 400 with 
the analytical estimate [ll| reported in Fig. [SJ; confirms the reliability of the numerical code. 

The stability of CPOs arising when a > a c — 8.34(1) Q is not known. However, the above arguments concerning 
the presence of a single degenerate band still hold, as confirmed by the spectrum plotted in Fig.[5jD, which corresponds 
to g = 0.4, and a — 9. This adds to the first band that is again located just below zero, with the exception of the first 
exponent which is exactly equal to zero, as a result of the quasi-periodic dynamics of the single neurons (the second 
zero Lyapunov exponent has been discarded while taking the Poincare section). 

Finally, in Fig. [5jl we have plotted the Lyapunov exponents of the first band multiplied by N 2 . The tendency 
to overlap of the spectra obtained for N = 50, 100, 200, and 400 indicate that we are again in the presence of a 
1/N 2 scaling as for the splay states, although finite-size corrections are more relevant in this case. Accordingl y, w e 
conjecture that the 1/N 2 dependence is a general property, which depends on the shape of the velocity field (see |lll|). 
rather than on the structure of the solution itself. 



B. Diluted Network 



Since the collective solutions observed in the globally coupled network are characterized by many weakly stable 
directions, it is reasonable to expect that generic perturbations of the network dynamics yield a chaotic evolution. On 
the other hand, in the thermodynamic limit we expect the inhomegenities induced by a random dilution to vanish. In 
fact, we have already seen (in the previous section) that a macroscopically periodic dynamics is eventually recovered. 

First of all we have verified that the network dynamics is chaotic for finite N, as soon as / > 0. In Fig. [f^, one 
can appreciate the growth of the maximal (positive) Lyapunov exponent with the fraction of missing links / (in the 
quenched case). 



8 




FIG. 6: (a) Maximal Lyapunov exponent Ai as a function of the percentage / of broken links for TV = 800 and a — 9. (b) 
Lyapunov spectrum {A^} versus k/N for a = 3 with networks of size TV = 50 (filled black circles) and TV = 100 (open red 
squares) and 20% of broken links. The data has been obtained by following the evolution of the event driven map and of the 
associated linearized equations ruling the evolution of the Lyapunov vectors for 3 — 4 x 10 s consecutive spikes, after discarding 
a transient composed by 300, 000 — 600, 000 spikes. The data refers to quenched disorder. 
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FIG. 7: First band of the Lyapunov spectra {At,} versus k/N for g = 0.4 with 20% of broken links for quenched disorder with 
a — 3 (a) and 9 (b) and for annealed disorder with a = 3 (c) and 9 (d). For the quenched disorder only spectra corresponding 
to system size TV = 50 are shown but for three different random network configurations, while in the annealed case spectra 
for network sizes TV = 50 (filled black circles) and TV = 100 (open red squares) are reported. The data has been obtained by 
following the evolution of the system and of its linearized copies for time lags similar to those reportd in the caption of Fig. [6] 



For a = 3, analogously to the globally coupled case, the spectrum is composed of two bands and two isolated 
eigenvalues (see Fig.|Bb). On the one hand, the degeneracy among the exponents lying in the second band is lifted (as 
a consequence of the disordered network structure) and the second band acquires a finite width. On the other hand, 
the comparison between the spectrum obtained for TV = 50 and TV = 100, reveals that the band width shrinks with 
TV around the value A = —a. This result again confirms that the inhomogenities of the disordered network vanish in 
the thermodynamic limit. A similar scenario occurs also in the CPO regime, as well as for annealed disorder. 

Furthermore, in Fig. [7] (where we plot the first band of the Lyapunov spectrum for quenched and annealed disorder, 
in both dynamical phases), we see a variable number of positive exponents, indicating that finite disordered networks 
are typically chaotic. More precisely, panels a and b of Fig. [3 refer to quenched disorder. In both cases, we present 
the spectra resulting from three different realizations of the disorder. We see that sample-to-sample fluctuations 
are quite relevant in the second part of the band, while they affect less the largest exponents. Another qualitative 
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observation concerns the shape of the spectrum that seems to be smoother in the presence of CPOs. Unfortunately, 
it is almost impossible to perform a quantitative scaling analysis of the spectrum (given the need to average over 
different realizations and to consider yet larger network sizes). In the annealed case, Fig. [7fc/d, there is no need 
to average over different realizations of the disorder and the Lyapunov spectra turn out to be smooth. However a 
scaling analysis is still beyond our computer capabilities: in both cases we report the spectra obtained for N = 50 
and N — 100 which are far from exhibiting a clean scaling behavior (for this reason we do not even dare to formulate 
a conjecture). 
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FIG. 8: Maximal Lyapunov exponent Ai versus N: (a) annealed case; (b) quenched case. In the quenched case also a power-law 
fit is reported (red dashed line) with decay exponent /3 = 0.51 ± 0.01. Fhe data has been obtained by following the evolution 

in the tangent space of the event driven map for 10 8 10 9 consecutive spikes, after discarding a transient of 200, 000 spikes. 

Both figures refer to 20 % of dilution, with parameters a = 9, g — 0.4 and a — 1.3. 



Given the difficulties of dealing with the whole Lyapunov spectrum, we have concentrated our efforts on the 
maximum exponent Ai, since we are thereby able to study larger networks. We limited ourselves to studying the 
CPO regime. As shown in Fig. [51 Ai decreases with iV as a power law: in the annealed case, Ai « 1/N, while in the 
quenched case Ai ~ 1/yN. This is at variance with the Kuramoto model, where in the "equivalent" quenched case, a 
1 /N behaviour has been observed [l6| . Accordingly, although annealed networks are more chaotic than the quenched 
ones over the numerically accessible network sizes, we conclude that the opposite will be true for yet larger networks. 



V. CONCLUSIONS AND PERSPECTIVES 



Our numerical analysis suggests that in the thermodynamic limit, a random, uncorrelated network behaves like 
a homogeneous globally coupled system with a rescaled value of the coupling constant to account for the different 
fraction of active links. This is because each neuron receives a random series of spikes: when the number of neurons 
increases, the number of spikes per unit time increases as well, while the relative fluctuations decrease and all neurons 
are affected by increasingly similar forcing fields (for a discussion on the thermodynamic limit in neural networks see 
also [TH, [13] )■ While there are little doubts that the collective motion is independent of the presence of disorder, less 
compelling is the evidence that the same is true for the stability properties. This is because numerical simulations are 
computationally expensive and it is not possible (at least for us) to study the scaling behavior of the entire Lyapunov 
spectrum in disordered systems. This is doable in homogeneous networks, thanks also to a faster convergence to the 
thermodynamic limit, and we have indeed observed that the first band of the Lyapunov spectrum scales as 1/N 2 , 
exactly as in the splay state, a case that is analytically known [ll|. In disordered networks we have nevertheless been 
able to study the scaling behaviour of the maximal Lyapunov exponent, finding that it is positive and scales as 1/N 13 , 
with P close to 1/2 (resp. 1) in the quenched (resp. annealed) case. This means that deterministic chaos disappears 
in the thermodynamic limit. This observation is consistent with the fact that the collective motion is periodic: a 
periodically forced phase oscillator (such as a LIF neuron) cannot be chaotic. However, nothing seems to prevent the 
collective motion itself to be chaotic. Whether different connection topologies have to be invoked or yet unidentified 
constraints forbid this to occur, remains as an open question. 
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